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ABSTRACT 


This  report  presents  a  first-order  discussion  of  the  perturbations 


in  the  orbital  motion  of  a  satellite  about  an  oblate  primary.  The 


approach  is  non-traditional  in  that  the  secular  changes  in  the  semi 


major  axis,  eccentricity,  inclination,  and  longitude  of  the  ascending 


node  are  all  derived  directly  from  the  translational  and  rotational 


equations  of  motion.  The  exactly  solvable  problem  of  motion  in  the 


primary's  equatorial  plane  is  used  to  obtain  the  advance  of  the  argument 


of  periastron  and  the  consequent  change  in  the  time  of  periastron  passage 


In  addition,  a  simple  discussion  of  a  very  restricted  three-body  problem 


is  given  to  illustrate  the  effects  of  third-body  perturbations 
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INTRODUCTION 


A.  Motivation  and  Outline 

Celestial  mechanics  is  one  of  the  oldest  fields  of  study  in  physics 
and  astronomy.  With  the  advent  of  the  space  age  and  the  rapid  development 
of  the  electronic  digital  computer,  the  field  has  undergone  a  resurgence 
of  both  theoretical  development  and  practical  use.  Old  subjects  have 

% 

their  own  vocabulary,  their  own  mathematical  methods,  their  own  symbolism, 

t  their  own  assumed  knowledge  base,  etc.  Although  many  scientists  and 

engineers  lack  this  specialized  training,  there  is  still  the  need  for 
them  to  be  cognizant  of  the  principal  results  of  the  two-body  problem 
and  of  perturbation  theory.  For  most  artificial  satellites,  the  principal 
perturbation  is  that  due  to  the  oblateness  of  the  earth  (symbolized  by  a 
non-zero  value  for  J ^) .  In  this  report  I  derive  all  of  the  first-order 
secular  variations  due  to  in  a  non-traditional  (for  a  celestial 
mechanic)  fashion.  Thus,  one  can  see,  without  the  appurtenances  of 
celestial  mechanics,  both  the  physical  cause  of  the  perturbations  and 
their  effects.  Some  of  these  results  are  scattered  in  the  standard 
reference  books.  However,  the  treatment  of  the  exact  solution  for  motion 
in  the  earth's  equatorial  plane  is  new.  Finally,  as  the  satellite's 
distance  from  the  earth  increases,  the  perturbations  of  its  motion  due 

■ 

to  the  sun  and  the  moon  become  more  important.  Hence,  a  very  simple 
treatment  of  third-body  perturbations  is  also  included. 

« 

The  report  has  some  caveats  the  reader  needs  to  know.  First, 
because  almost  no  new  physics  is  involved,  the  equations  have  been  kept 
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to  a  minimum.  Second,  to  keep  th  framewor.  ot  die  Jiacussion  constant, 
this  is  a  geocentric  report.  Nowhere  is  any  particular  property  of  the 
earth  used  nor  does  the  geocentric  aspect  limit  the  discussion.  Lastly, 
except  for  Si  ID*,  the  entire  report  is  for  first-order  perturbation 
theory  only.  Extrapolating  the  conclusions  reached  herein  beyond  this 
limitation  is  not  necessarily  correct. 

We  next  present  brief  summaries  of  two-body  elliptic  motion  and  the 
geopotential.  Following  that  is  a  detailed  discussion  of  motion  in  the 
earth's  equatorial  plane.  This  illustrates  various  techniques  of  analysis 
and  an  analytic  solution  to  such  a  problem.  We  then  turn  to  motion  in 
any  plane  and  exhibit  another  method  of  analysis.  This  time  the  equations 
for  the  rates  of  change  of  the  angular  momentum  and  the  energy  are  used 
to  further  our  understanding  of  the  motion.  This  section  closes  with 
the  perturbation  equations  in  the  classical  form.  The  last  section 
deals  with  a  primitive  three-body  problem  so  that  the  effects  of  the  sun 
and  moon  can  be  illustrated. 

B.  Brief  Summary  of  Two-Body  Elliptic  Motion 

The  classical  two-body  problem  has  the  artifical  satellite,  of 
negligible  mass  ra,  revolving  about  the  earth  whose  mass  is  M  .  The 
coordinate  system  with  origin  at  the  center  of  mass  of  the  earth  is 
assumed  to  be  an  inertial  one.  We  will  take  the  extension  of  the  earth’s 
equator  as  the  xy  plane  and  have  the  x  axis  point  in  the  direction  of  the  Vernal 
Equinox.  The  z  axis  points  in  the  direction  of  the  North  Celestial  Pole 
♦Equations  (46,  47)  are  exact  but  our  handling  of  them  is  first-order. 


and  the  y  axis  completes  a  right-handed  triple.  The  force  is  central 
[the  potential  is  given  in  Eq.  (11c)]  so  the  motion  takes  place  in  a 
plane.  The  plane's  normal  is  the  angular  momentum  vector.  The  form  of 
the  orbit  is  an  ellipse  which  can  be  parametrically  described  by 
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r  =  a(l-e  )/(l  +  ecosv)  (la) 


or 


r  *  a(l  -  ecosE) .  (lb) 

The  ellipse's  semi-major  axis  is  a,  its  eccentricity  is  e,  v  is  the 
azimuthal  coordinate  in  the  orbital  plane  measured  from  the  point  of 
closest  approach  (i.e.,  perigee;  v  is  also  called  the  true  anomaly),  and 
E  is  an  auxiliary  angular  variable  named  the  eccentric  anomaly.  For  an 
arbitrary  orientation  of  polar  coordinates  in  the  orbital  plane  the 
distance,  r,  would  be  given  by 

r  ■  a(l-e2)/[l  +  ecos(u-w)],  (lc) 

u  ■  v  +  w,  (2) 

where  u>  is  called  the  argument  of  perigee.  See  Figure  1.  The  time 
dependence  of  the  orbit  is  Implicitly  given  by  the  relationships  between 
the  true  anomaly,  the  eccentric  anomaly,  and  the  mean  anomaly,  M; 
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. _ 


(3) 


tan(v/2)  =  1(1  +  e)/(l  -  e) "tan(E/2) , 

E  -  esinE  -  M  =  n(t  -  T)  -  n(t  -  T  )  +  M  ,  (4) 

o  o 

where  the  mean  motion,  n,  is  related  to  the  semi-major  axis  and  period, 

P,  of  the  motion  by 

.  1/2  -  1/2 

n  =  (GM^/aJ)  =  (y/aJ)  =  2u/P.  (5) 

G  is  the  universal  constant  of  gravitation.  The  standard  epoch  is  the  time 

of  perigee  passage,  T,  when  M  =  0.  If  some  other,  arbitrary,  epoch  t  =  Tq 

is  used,  then  M  =  M  at  that  instant, 
o 

To  complete  the  three  dimensional  picture  of  the  motion,  we  need  to 

A 

know  the  direction  of  the  angular  momentum  vector,  L, 

A 

L  »  (sinftsini,  -cosfisini,  cosi) .  (6a) 

The  inclination  of  the  orbital  plane  to  the  equator  is  symbolized  by 
ie[0,180°].  That  place  on  the  equator  where  the  plane  of  the  orbit 
intersects  it  and  the  satellite  is  moving  northward  is  called  the  ascending 
node.  Q  is  the  longitude  of  the  ascending  node,  fie[0,360°).  See  Figure  2. 
The  magnitude  of  the  angular  momentum  vector  is  given  by  mL, 

j  1/2 

L  =  [pa(l  -  e  )]  ,  L  *  LL.  (6b) 
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The  angle  (for  which  there  is  no  standard  notation)  w  +  (2  is  called  the 
longitude  of  perigee.  Finally,  to  go  from  the  cartesian  orbital  plane 
location  vector 


_r  =  r(cosv,  sinv.o)  , 


(7a) 


to  the  cartesian  geocentric  location  vector  or  spherical  geocentric 
location  vector. 


r  = 


/ cos6cosa  \ 
r|  cosdsina 
\  sin6  / 


(7b) 


we  need  to  rotate  by  the  orthogonal  matrix  S, 
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(7c) 


The  results  are 


x/r  ■  cosHcosu  -  sinftcosisinu. 


y/r  -  sinflcosu  +  cosftcosisinu 
z/r  ■  sinisinu. 


(8) 


or 


a  =  9.  +  tan  (cositanu) , 


(9a) 


<5  =  sin  (sinisinu). 


(9b) 


where  u  was  defined  in  Eq .  (2).  Alternative  expressions  of  the  invariable 
plane  are 


xsinftsini  -  ycosftsini  +  zcosi  =  0 


(10a) 


and 


tan6  =  tanisin(a  -  £i) .  (10b) 

All  of  this  is  derived  directly  from  the  equations  of  motion.  In  the 
plane  of  the  orbit  their  form  is 


r  -  rv2  *  -  3U/3r;  dw/dt  =  w  Vw,  (11a) 

rv  +  2rv  -  -  (l/r)3U/3v,  (lib) 


U  -  -y/r. 


(11c) 


The  key  to  the  solution  of  the  two-body  problem  and  the  remainder  of  the 
analysis  presented  herein  is  the  ability  to  analytically  integrate  Eq.  (lib) 
once. 
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C.  Brief  Summary  of  the  Geopotential  Through  J2 
The  gravitational  potential  at  a  point  located  r_  outside  of  a  body 
whose  mass  density  is  p  is  given  by 


U(r) 


-G  J 
Volume 


p(r')dr' 


r  -  r' 


(12) 


If  one  uses  the  spherical  harmonic  expansion  of  the  denominator  of  the 
integral  one  finds,  through  second-order  terms. 


U(r)  =  (G/r)/p(r')dr’  -  (G/r2)/r 'p(r ')cosgdr ’ 


+  (G/2r2)/(r')2p(r.')(l  -  3cos2g)dr\ 


(13) 


The  angle  between  _r  and  r'  is  denoted  by  g.  The  first  integral  is  the  total 
mass  (M^  for  the  earth) .  The  second  Integral  vanishes  because  the  origin  of 
the  coordinate  system  is  at  the  center  of  mass.  The  third  integral  is  equal 
to  [(three  times  the  moment  of  inertia  about  the  line  joining  the  center  of 
mass  and  r )  minus  (the  sum  of  the  principal  moments  of  inertia)].  If  we  employ 
the  symmetry  of  an  oblate  spheroid  and  express  the  last  integral  in  terras  of 
the  principal  moments  of  Inertia  of  the  spheroid  (A  -  equatorial  moment  of 
inertia,  C  -  polar  moment  of  inertia)  then 

U(r,6)  -  -  (y/r) [1  +  J2R2(1  -  3sin26) / (2r2) ]  (14a) 

where 

J2  -  (A  -  C)/(M^R2)  -  1.0827  x  10~3  for  the  earth,  (14b) 
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and  6  is  the  geocentric  latitude  of  i-.  We  shai.V  abbreviate 

n  -  j2R9>  U5) 

and  note  that  the  ratio  of  the  forces  due  to  the  earth's  total  mass  and  that 
due  to  its  oblateness  on  an  object  a  distance  R  >  R^  away  is 

e  *  n/R2.  (16) 

For  a  stationary  satellite,  e  “  2.5  x  10  ^  -  5  arc  seconds. 
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II.  MOTION  IN  THE  EARTH'S  EQUATORIAL  PLANE 


With  6  »  0  the  potential,  from  Eq.  (14a),  is 

U(r,0)  -  -  (p/r)(l  +  J2R^/2r2)  -  -  (p/r) (1  +  n/2r2) .  (17) 

The  equations  of  motion,  Eqs.  (11a,  lib),  have  the  two  integrals 

(r2  +  r2v2)/2  +  U(r,0)  =  E(<  0  here),  (18a) 

2 , 

r  v  *  L  (>_  0  always).  (18b) 

The  existence  of  the  angular  momentum  integral  allows  a  discussion  of  the 
orbit,  r(v),  instead  of  the  time  dependence  of  the  motion.  Since, 

r'  =  dr/dv  -  r/v,  r  -  Lr'/r2.  (19) 

The  conservation  of  energy  equation,  Eq.  (18a),  can  be  written  as 

(r’)2  +  r2  -  2r4[E  -  U(r,0)]/L2.  (20) 

When  J2  ■  0  the  semi-major  axis  and  eccentricity  were  defined  by 


As  long  as  is  small,  those  quantities  will  have  a  meaningful  interpretation 

as  a  "semi-major  axis"  and  "eccentricity".  In  any  case,  they  too  are  constants 

2 

of  the  motion.  Eliminating  E  and  L  from  Eq.  (20),  the  equation  describing 
the  orbit  may  be  reexpressed  as 

(r1)2  -  -  rF(r)/[a2(l  -  e2)]I/2,  (22a) 

with 


F(r)  «  r (r  -  r+)(r  -  r_)  -  an. 


(22b) 


and 


r  ■  a(l  +  e) .  (22c) 

This  will  be  the  starting  point  for  the  remainder  of  the  analysis  presented 
in  8 II.  We  note  that  if  Jj  were  zero,  then  the  roots  of  F  would  be  0,  r_ 
(perigee  distance),  and  r+  (apogee  distance). 

A.  Circular  Orbits 

A  circular  orbit  is  always  possible  for  non-zero  an<*  always  stable 

for  small  J ^ •  Since  r'  “0,  the  radius  of  the  orbit  is  a  root  of  F.  It's 

2 

actually  simpler  to  derive  it  from  Eq.  (20).  The  result  is  (e  ■  r|/a  ) 

radius  -  (a/2) [1  +  (1  -  6e)1/2]<a.  (23) 
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Thus,  the  only  effect  is  to  reduce  the  size  of  the  orbit.  But  from  Eq.  (17)  we 


see  that  if  r  is  constant,  then  a  positive  J ^  merely  augments  the  central  mass; 
hence  for  the  same  energy  the  satellite  would  be  more  tightly  bound. 

B.  "Elliptical  Orbits":  Apogee  and  Perigee 

It  is  clear  from  Eq.  (22a)  that  the  roots  of  F  represent  the  turning 

points  of  the  orbit.  It  is  a  simple  matter  to  show  that  F  has  three  real 

roots  (two  of  which  are  equal  for  a  circular  orbit)  and  we  know  that  as 

^2  0  they  must  approach  0,  r  ,  and  r+.  The  one  near  zero  is  not  a  physical 

turning  point.  Let  us  call  the  three  true  roots  of  F  (for  J_  4  0) ,  R  , 

i  o 

R  (for  perigee),  and  R  (for  apogee).  Their  values  can  be  approximated 
P  « 

2 

through  the  use  of  Newton's  method.  One  finds  (e  =  n /a  still) 


Rq  =  ae/(l  -  e  ), 


(24a) 


R  =  r_  -  ae/[2e(l  -  e)]  <  r_, 


(24b) 


Rfl  =  r+  +  ae/I2e(l  +  e)]  >  r+. 


(24c) 


Thus,  perigee  is  closer  in  and  apogee  further  out  for  the  perturbed  "ellipse" 
than  for  the  unperturbed  ellipse.  The  "semi-major  axis".  A,  and  "eccentricity", 
E,  can  be  derived  by  using  the  formulas  for  elliptical  motion,  namely 


2A  -  R  +  R  , 
a  p* 


(25a) 


■» 


(25b) 


2AE  *  R  -  R  , 
a  p’ 


or 


A/a  -  1  -  e / [ 2 (1  -  e2)]  <  1,  (26a) 

E/e  *  1  +  e/(l  +  e2)/[2e2(l  -  e2)]  >  1.  (26b) 

We  already  knew  that  the  "focal  distance"  was  smaller  and  we  now  see  that 
the  reason  is  that  the  "eccentricity"  is  larger  and  the  "semi-major  axis"  is 
smaller. 

It  now  seems  logical  to  ask  about  the  utility  of 

r  -  A(1  -  E2) / (1  +  Ecosv),  (27) 

as  a  representation  of  the  orbit.  If  we  express  this  through  first  order  in 
H  we  find 


2  3 

,  »  ae(cosv  +  3e  +  3e  cosv  +  e  ) 
r  -  r  (v)  -  — 1 - 5 - o  * 

2e(l  -  e  )(1  +  ecosv) 


where 


r  (v)  ■  a(l  -  e  )/(l  +  ecosv). 
o 


(28a) 


(28b) 
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Equation  (27)  can  only  represent  the  orbit  except  for  those  terms  in  v  which 
vanish  at  v  -  0  and  v  ■  tt.  However,  the  value  of  r  at  the  turning  points  is 


independent  of  v,  so  that  there  can  be  no  secular  change  in  A  or  E  due  to  J 


Elliptical  Orbits":  The  Orbit 


Assuming  that  Eqs. (22)  can  not  be  exactly  solved,  we  try  a  physicist's 
perturbation  approach.  We  write 


substitute  this  into  the  equation,  separate  powers  of  e,  and  solve  the  result 


ing  system  of  equations.  If  we  do  this  [but  we  actually  use  r(v)  *  r  (v) 

o 

2 

+  e(l  -  e  )r.(v)  for  simplicity]  the  coefficient  of  e  (which  must  vanish)  is 


P(v)  ■  (ecos  v  -  cosv  -  2e)cosv/(l  +  ecosv) 


q (v)  ■  a(l  +  ecosv) cscv/[2e(l  -  e  ;  ] 


The  solution  of  Eqs.  (30)  for  r. (v)  is 


Csinv 


alcosv  +  3e  +  3e 


cosv  +  vsinv)  +  e  (1  +  sin  v 


(1  +  ecosv) 


2e(l  -  e  )(1  +  ecosv) 


The  arbitrary  constant  of  integration  is  C  and  the  initial  condition 


is  r(0)  •  R  with  the  expression  (24b)  used  for  R  .  Although  the  initial 
P  P 

condition  does  not  determine  the  value  of  C,  the  existence  and  uniqueness 
theorem  for  the  differential  equation  (30a-c)  allows  us  to  set  C  *  0.  Hence, 


r  (v) 


r  (v)  - 
o 


ae[cosv  +  3e  +  3e  (cosv  +  vsinv)  + 


2e(l  - 


e2)(l 


e~*(l  +  sin2v)  ] 


(31) 


+  ecosv)' 


Compare  with  Eqs.  (28). 

We  can  now  extend  our  analysis  to  compute  the  advance  of  perigee. 

The  presence  of  the  secular  term,  oc  vsinv,  means  that  although  r(0)  “  t(2tt) 

■  r(4iT)  =  .  .  .  and  r(n)  ■  r(3u)  ■  r(5u)  .  .  .  (cf.  page  14),  it  is  not 
necessarily  true  that  the  minimum  and  maximum  values  of  r  occur  at  these  true 
anomalies.  By  construction,  the  first  perigee  occurs  for  v  «  0.  We  know  the 
first  apogee  must  occur  near  v  =  it  [cf.  Eq.  (la)  or  (28b)].  By  computing 
r' (v)  from  Eq.  (31),  setting  v  »  v  -  tt  +  eV  in  the  result,  and  solving  for 

3  3 

by  expanding  all  terms  through  first  order  in  e  we  find 


vfl  *  tt { 1  +  3e/[2(l  -  e2)2]}  >  it. 


(32a) 


Similarly,  the  second  perigee  must  occur  near  v  «  2ir,  and  an  analogous 
procedure  yields 


v  -  2v  >  2-n .  (32b) 

P  a 


Hence,  the  advance  of  the  argument  of  perigee  is  symmetrically  distributed 
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over  each  half  of  the  orbit.  Clearly  the  time  of  perigee  passage  does  not 


increase  by  only  P.  Its  change  is  given  by  (vp  -  2tt)/<v>  where  the  angular 
bracket  denotes  an  average  angular  rate.  A  simple  average  for  v  is  obtained 
from  Eq.  (18b)  at  the  perigee  and  apogee  points,  whence 

T  -  P  ®  3ire/[n(l  -  e2)3/2]  >  0.  (33) 

The  reader  may  notice  that  r(ir)  =  R  with  R  given  by  Eq.  (24c)  and  then 

3L  cl 

wonder  why  v  of  Eq.  (32a)  is  greater  than  ir.  The  reason  is  that  the  Taylor 
series  for  the  cosine  contains  no  first-order  terms. 

Finally,  to  illustrate  these  effects  in  an  exaggerated  fashion  Eq.  (31) 
has  been  graphed  for  a  *  1,  e  =  0.25,  e  =  0.05  and  e  =  0.10  at  10°  intervals 
in  the  true  anomaly.  The  solution  of  Eqs.  (22)  in  the  form  given  by  Eq.  (29) 
can  not  be  reliably  carried  out  beyond  ev  ^  1  except  by  piecewise  continuation. 
D.  The  Exact  Solution 

The  exact  solution  of  Eqs.  (22a,  b)  involves  elliptic  functions.  We 
need,  at  minimum,  the  sine  amplitude,  sn,  the  cosine  amplitude,  cn,  and  the 
delta  amplitude,  dn.  We  can  define  the  sine  amplitude  via 

p  ■  /  (1  -  K2sin20)^2d9,  sn(p,  k)  =  sin<(>,  (34b) 

o 

and  the  others  by 

sn2(P,  <)  +  cn2(p,  K)  -  1, 

2  2  9 

k  sn  (p,  tc)  +  dnz (p ,  tc)  .  i. 
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(34a) 

(34c) 


lUI-S-tl 


Fig.  4.  Perturbed  motion  (solid  curve)  and  unperturbed  motion  (dashed 
curve)  for  a  =  1,  e  =  0.25,  e  =  0.05  and  a  =  1,  e  =  0.25,  e  =  0. 


While  these  functions  have  various  addition  theorems,  half-angle  formulas, 
periodicities,  etc.  we  shall  need  few  of  them  here.  In  particular  though. 


sn2(y/2,<)  -  U  -  cn(y,<)]/[l  +  dn(y,<)], 

2 

cn  (y/2,<)  -  Icn(y.ic)  +  dn(y,ic)]/[l  +  dn(y,ic)), 

dn2 (y/2 ,<)  *  Icn(y,<)  +  dn(u,<)]/Il  +  cn(y,K)], 

sn(u  +  v.k)  -  Isn-(y><)cn<^.<)dn(v.<)  +  sn(v,<)cn(y ,K)dn(y,K)] 

2  2  2  * 
1  -  k  sn  (y,K)sn  (v,k) 

sn[K(<) ,k]  *  1, 

sn[y  +  4K(k),  k]  =  sn(y,K), 


(35a) 

(35b) 

(35c) 

(35d) 

(36a) 

(36b) 


where  the  complete  elliptic  integral  of  the  first  kind,  K,  is  given  by 

K(k)  -  /  /  (1  -  K2sin20)1/2d6.  (36c) 

o 

The  quantity  K  is  called  the  modulus.  The  last  thing  we  need  to  know  is  that 
any  integral  of  the  form 

R  (x)  +  R  (x)P1/2(x) 

1  "  /  - wo -  d*.  (37) 

yx  R3(x)  +  RA(x)PA/z(x) 

where  R^(x),  i  ”  1,2, 3, 4  is  a  rational  function  of  x  and  P  is  a  quartic 
polynomial  in  x  with  real  coefficients,  can  be  reduced  to  an  elliptic 
function  (or  sum  of  them) . 
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Now,  Eqs.  (22a,  22b)  are  of  the  form  (37)  if  we  rewrite  them  as 


dr/dv  -  +  {-rF(r)/Ia2(l  -  e2)]}1/2. 


(38) 


We  also  know  that  r'  >  0  for  ve[0,  v  ]  and  r'  <  0  for  ve[v  ,  v  ].  The  true 

—  *  a  —  a  p 

anomalies  of  first  apogee  passage  and  second  perigee  passages  are  approximately 
given  in  Eqs.  (32).  Thus, 


ds 


+1 


Rp  [-sF(s) ] 


1/2 


a(l  -  e)  o 


/  dw;  R  £r£R,0<.v£v,  (39a) 


2^1/2  -p  -  -  -  -  -  -  a' 


ds 


-1 


.1/2 


R  [-sF(s) ]  a(l  -  e)  va 

oi  a 


/  dw;  R  <r<R,v  <v<v.  (39b) 

'  9  n  —  —  a  '  a  —  —  n 


2n1/2  '  -p  -  -  -  &’  a  -  -  p 


There  exists  a  standard  technique  for  reducing  integrals  of  the  form  (37) 
to  the  standard  forms  [of  which  (34a)  is  one  example] .  The  results  are 
for  our  cases 

1/2 


yv/2  *  sn 


•< 


R  -  R  r  -  R 

a  o  .  p 


,R  -  R  r  -  R 

Lap  o 


Y(v  -  va)/2  ”  an" 


-(Li. 

\L a  p 


R  -  r 
a 


4 
,  ; 


ve[0,va] , 


ve[va,vp]. 


(40a) 


(40b) 


where  the  modulus,  k,  is  given  by 
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and  the  quantity  y  ls  given  by 


transformed  into 


R  R  {[1  +  dn(yv.k)]  -lc  [1  -  cn(yv,k)]} 


We  can  now  rederive  all  of  our  earlier  work  as  a  check.  From  Eq.  (36c) 


Then  from  Eq.  (36a)  it  follows  that 


it  {1  +  3e/[2(l  -  eVl) 


From  Eqs.  (39)  and  this  result  Eq.  (32b)  follows  immediately.  Also,  from 


Eqs.  (34)  it  follows  that  if  K  is  small 


2 

cn(y,<)  -  cosy  +  (k  /4)(y  -  sinycosy)siny 


(34d) 


dn(y,K)  -  1  -  (<2/2)sin2y. 


(34e) 


so,  from  these  approximations,  Eqs.  (40c,  d  and  41)  the  result  in  Eq.  (31) 
can  be  derived.  The  periodicity  of  the  motion  is  assured  because  cn  and 
dn  share  4K  as  a  period.  The  only  remaining  point  is  the  time  dependence 
of  r  and  v  which  can  be  computed  directly  from  Eqs.  (18)  and  involves 
other  elliptic  functions. 


111.  GENERAL  MOTION 


In  this  section  the  rates  of  change  of  the  energy  and  the  angular 
momentum  are  used  to  derive  the  first-order  secular  perturbations  in 
four  (a,  e,  i,  and  ft)  of  the  six  orbital  elements.  The  analysis  of  the 
preceding  section  has  determined  the  changes  in  the  other  two  except  for 
a  multiplicative  dependence  on  inclination  (which  is  non-zero  for  i  -  0)  and 
an  additive  dependence  on  inclination  (which  vanishes  for  i  -  0) .  At  the 
end  of  this  section  one  version  of  the  classical  perturbation  equations  is 
given  and  solved. 

The  potential  was  given  in  Eqs.(14,  15). 

A.  The  Work  Done  by  3^ 

The  work  done  by  Jj,  dE/dt,  is  given  by 

dE/dt  -  -  VU  ‘  r.  (42a) 

Over  the  course  of  one  revolution  the  change  in  E  is 

A E  -  §Edt .  (42b) 

One  can't  simply  note  that  V  x  VU  «  0  and  use  Green's  theorem  to  evaluate 
this  though.  The  sign,  §,  means  integrate  around  one  complete  circuit  of 
the  perturbed  orbit.  We  can  not  assume  that  this  leads  to  the  same  point 
in  space.  If  we  split  U  into  Uq  +  r|U,  where  U  ■  -  p/r  then. 
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AE  «  AE  +  AE.  -  §[  -  VU  -  nVU.]  *  rdt. 
o  1  o  1  — 


(42c) 


Now,  to  first-order  in  n  we  can  deform  the  contour  for  the  Uq  integration 
to  be  the  unperturbed  orbit.  Now  we  exploit  the  fact  that  the  force  due  to  UQ 
is  conservative  so  AEq  ■  0.  Since  f|U^  is  already  first-order,  we  can  again 
deform  the  contour  and  AE^  »  0.  But, 


E  -  -  y/2a, 


(21) 


so 


Aa  -  0  or  da/dt)  *  0,  (43a) 

sec 

where  da/dt)8gc  is  the  first-order  secular  rate  of  change  of  the  semi-major 
axis  due  to  J2> 

B.  The  Torque  Due  to  Jj 

The  rate  of  change  of  the  angular  momentum  due  to  is  given  by 

dL/dt  “txF«-rxVU»  -  nrxVU^,  (44a) 

since  Uq  represents  a  central  force.  By  direct  computation  we  find 

L  -  -  3(ynz/r5)(y,  -x,  0),  (44b) 

and 
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AL  -  /  Ldt  -  -  (cosfl,  sinft,  0) 

o  2naZ(l  -  e  )J/ 

-  (ALx,  ALy,  ALz). 

From  Eqs.  (6) 

AL  -  (L  +  AL)(sinln  +  Aft]sin[i  +  Ai], 

-  cos[ft  +  Aft]sin[i  +  Ai],  cos[i  +  Ai])  -  L 

-  AL(sinftsini,  -  cosfisini,  cosi) 

+  LAiHcosftsini,  sinftsini,  0) 

+  LAi(sinftcosi,  -  cosftcosi,  -  sini) . 

From  Eqs.  (45)  we  find 

AL  ■  AL  sln&sinl  -  AL  cosftsini  +  AL  cosi, 
x  y  z 

LAflsini  ■  AL  cosfl  +  AL  sinft, 
x  y 

LAi  *  AL  sinftcosi  -  AL  cosftcosi  -  AL  sini, 
x  y  z  * 

so,  after  dividing  by  P  and  setting  Ai/P  ■  di/dt)  et 


dQ/dt) 


-  3encosi/[2(l  -  e2)2] 


(43c) 


sec 

da/dt)  -  [2ea/(l  -  e2)]de/dt)  ,  (46) 

see  sec 

where  Eq.  (6b)  has  been  used  to  derive  Eq.  (46)  from  dL/dt)  =  0.  Equation 

sec 

(46)  and  Eq.  (43a)  lead  to 

de/dt)sec  =  °‘  (43d) 

Equations  (43)  then  summarize  our  results.  We  see,  therefore,  that  the 
application  of  the  equations  for  the  time  rates  of  change  and  E  and  L  have 
simply  and  directly  enabled  us  to  solve  two-thirds  of  the  perturbation 
problem.  This  approach  also  gives  a  direct  physical  insight  into  the  cause 
of  the  perturbations  as  well  as  their  effects. 

C.  Classical  Perturbation  Theory 

The  two-body  problem  with  the  potential  given  in  Eq.  (11c)  forms  the 
cornerstone  of  celestial  mechanics.  It  represents  the  situation  exactly  and 
can  be  solved  exactly.  Hence,  it  is  treated  in  all  mechanics  texts  and  is 
fully  developed  analytically.  The  coordinate  system  represented  by  the  six 
orbital  elements  has  become  familiar  and  useful.  Intuitively,  if  the  total 
real  force  in  a  problem  only  departs  slightly  from  that  given  by  -VUo,  the 
orbital  elements  themselves  can  only  change  slowly.  This  is  the  foundation 
of  classical  perturbation  theory.  The  three,  second-order,  differential 
equations  given  by  F  »  rna  are  replaced  by  six,  first-order  differential 
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equations  for  de/dt,  di/dt,  etc.  This  can  be  performed  In  a  variety  of 
ways  for  both  non-conservative  and  conservative  forces.  The  reduction 
is  purely  a  computational  problem  and  can  be  simplified  by  using  the 
Lagrangian  formalism.  The  important  restriction  is  the  condition  of 
osculation.  This  mathematical  statement's  content  is  that  both  £  and  £ 
can  be  computed  from  the  osculating  values  of  the  orbital  elements  via 
the  ordinary,  two-body,  Keplerian  formulas. 

I  shall  give  one  form  of  the  conservative  equations  here.  Others, 
to  be  found  in  the  standard  reference  texts,  will  vary  in  the  choice  of 
orbital  parameters.  Let  the  perturbing  force  be  given  by 


F 

—pert 


-  +  VU' 


(46) 


Then, 


da/dt  -  [2/(na)J3U'/3M,  (47a) 

de/dt  -  {(1  -  e2)/[na2e]}3UV3M  -  {(1  -  e2)1/2/[na2e]}3U73u>,  (47b) 

dw/dt  -  -  {coti/[na2(l  -  e2)1/2]}3U’/3i  +  {(1  -  e2)1/2/[na2e]}3U’ /3e,  (47c) 
di/dt  -  {coti/[na2(l  -  e2)l/2]}3U' /3u>  -  {csci/[na2(l  -  e2)1/2]}3U' /3ft, (47d) 
dft/dt  -  (cscl/[na2(l  -  e2)1/2]}3U'/3i,  (47e) 
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dM/dt  -  -  {(1  -  e2)2/[na2e]}3U’/3e  -  I2/(na)]3U'/3a.  (47f) 

As  an  illustration  of  how  to  use  these,  we  take 

U’  -  (ny)(l  -  3sin26)/(2r3),  (48) 

from  Eq.  (14a).  Since  U'  is  already  of  first  order,  we  can  replace  6  and  r 
via  Eqs.  (9b  and  1)  in  U',  compute  the  derivatives  on  the  right-hand  sides  of 
Eqs.  (47),  regard  a,  e,  ...»  as  constants,  and  integrate.  This  will  give  us 
knowledge  we  have  not  yet  derived,  the  explicit  time  dependence  of  a,  e,  .... 
as  opposed  to  their  secular  rates  of  change.  Of  course,  if  we  integrate  over 
a  period  we  reproduce  Eqs.  (43a  -  d)  and  can  derive  the  inclination  dependence 
for  du)/dt  and  dT/dt.  The  latter  results  are 

<Wdt)gec  -  3en(5co82i  -  1)/I4(1  -  e2)2],  (43e) 

dT/dt) sec  -  -  3e(l  -  3sin2isin2w)/[2(l  -  e)3].  (43f) 

Note  that  for  cos2i  -  1/5  (i  -  63°  26'  5'.'82  or  i  -  116°  33'  54’.'18)  there 
is  no  perigee  advance  due  to  (in  first-order  theory!). 

Equations  (47),  or  their  non-conservative  generalization,  give  us  a 
very  powerful  tool  for  studying  long-term  perturbing  effects  analytically 
since  short-period  perturbations  can  be  left  out  (they'll  be  averaged  over 
anyway) . 
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IV.  A  VERY  RESTRICTED  THREE-BODY  PROBLEM 


We  could,  of  course,  discuss  third  body  perturbations  within  the  frame¬ 
work  of  classical  perturbation  theory.  However,  in  its  simplest  formulation, 
this  problem  admits  an  analytical  first  order  solution  so  we'll  follow  a 
physicist's  approach.  Moreover,  this  time  we  will  solve  directly  for  the 
time  dependence  of  r  and  v  instead  of  for  r(v).  To  set  the  stage  the 
satellite's  unperturbed  orbit  is  circular,  viz 

r  -  a,  v  =  2irt/P  with  Ai^a^/P2  »  y.  (49a) 


In  the  plane  of  the  orbit  lies  another  body,  of  mass  Mp, 
about  the  earth  in  a  circular  orbit. 


also  revolving 


rp  ‘  V  VP '  2n,\  *  *•  p '  “•  <«">> 

However,  a^»a  or  P»P^.  Thus,  we  can  regard  the  perturbing  body  to  be  fixed 
for  a  few  revolutions  of  the  satellite  (a  poor  approximation  for  the  Moon  but 
a  reasonable  one  for  the  Sun).  The  potential  for  the  problem  is 


U  -  -y/r  -  GM/ | r  -  r|, 
P  P 


(50) 


or 


U 


-y/r  -  GM  /la2  +  r2  -  2ra  cos(v  -  v  J1^2. 
P  P  P  P 


(51) 
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The  equations  of  motion  of  the  satellite  are  still  Eqs.  (11a,  b) .  If  we 
use  the  approximations  mentioned  above  they  become 


2  2 

r  -  rv  «  -y/r  +  yecos(v  -  if>) , 

rv  +  2rv  -  -  esin(v  -  <J/) , 

e  -  (GM/a2) /GM/a2) . 

P  P  * 

For  the  Sun  e  *  0.026,  for  the  Moon  e  ■  0.00015.  We  set 
r  -  a  +  er^,  v  «  27rt/P  +  ev.^, 
substitute  into  Eqs.  (52),  linearize,  and  find 

-  12ir2r^/P^  -  4?rav^/P  ■  (y/a2)cos(2irt/P  -  i|»), 
av^  +  Uvifl  -  -  (y/a2)sin(2irt/P  -  \|i) . 

As  initial  conditions  we  use 

r(0)  -  0,  v(0)  -  0, 

*(0)  -  0,  *(0)  -  2n/P. 


(52a) 

(52b) 

(52c) 

(53) 

(54a) 

(54b) 

(55a) 

(55b) 
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We  can  integrate  the  equation  once,  substitute  that  into  the 
equation  and  then  integrate  it  twice.  We  then  return  to  the  equation 
to  find  it.  The  answer  is 

r/a  ■  1  +  e[Bcos(2irt/P  +  $)  +  (3irt/P)sin(2irt/P  -  40  -  2cos40.  (56a) 

v  »  2nt/P  -  e{2Bsin(2irt/P  +  40  -  (3itt/P)  [cos(2irt/P  -  40  +  cos4>] 

+  2sin(2irt/P  -  ^)  +  5sinij;},  (56b) 


where 


Bcosi{)  -  2cos<|), 


(56c) 


Bsin$  ■  -(3/2)8in\|<, 


(56d) 


and  this  solution  is  valid  for  times  |t|  <  P/(2itc). 

Given  the  artificial  nature  of  the  problem,  it  seems  best  just  to 

illustrate  the  motion  for  a  few  values  of  e  and  <p.  Figures  5-8  show  the 

unperturbed  circular  orbit  as  a  full  curve  and  the  perturbed  orbit  as  a 

series  of  equally  spaced  dots  [At  ■  (1/36) 'th  of  the  unperturbed  period]  for 

i|/  ■  0(90)360°,  e  "  0.05  from  t  ■  0  until  t  ■  P.  From  the  location  of  M  and 

P 

the  direction  of  the  Initial  velocity  vector,  we  can  simply  interpret  the 
results.  It  is  a  numerical  accident  that  for  4*  ■  90°  or  270°  v(P)  ■  0. 
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This  report  presents  a  first-order  discussion  of  the  perturbations  in  the  orbital  motion  of  a 
satellite  about  an  oblate  primary.  The  approach  Is  non-tradltlonal  In  that  the  secular  changes  In 
the  semi-major  axis,  eccentricity,  Inclination,  and  longitude  of  the  ascending  node  are  all  derived 
directly  from  the  translational  and  rotational  equations  of  motion.  The  exactly  solvable  problem 
of  motion  In  the  primary's  equatorial  plane  is  used  to  obtain  the  advance  of  the  argument  of  perlastron 
and  the  consequent  change  In  the  time  of  perlastron  passage.  In  addition,  a  simple  discussion  of  a 
very  restricted  three-body  problem  is  given  to  illustrate  the  effects  of  third-body  perturbations. 
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